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We examine the effects of disorder on dimerized quantum antiferromagnets in a magnetic field, 
using the mapping to a lattice gas of hard-core bosons with finite-range interactions. Combining a 
strong-coupling expansion, the replica method, and a one-loop renormalization group analysis, we 
investigate the nature of the glass phases formed. We find that away from the tips of the Mott lobes, 
the transition is from a Mott insulator to a compressible Bose glass, however the compressibility at 
the tips is strongly suppressed. We identify this finding with the presence of a rare Mott glass phase 
not previously described by any analytic theory for this model and demonstrate that the inclusion 
of replica symmetry breaking is vital to correctly describe the glassy phases. This result suggests 
that the formation of Bose and Mott glass phases is not simply a weak localization phenomenon but 
is indicative of much richer physics. We discuss our results in the context of both ultracold atomic 
gases and spin-dimer materials. 

PACS numbers: 05.30.Jp, 64.70.P-, 64.60.ae, 64.70.Tg 

The disordered Bose-Hubbard model is an ideal sys¬ 
tem for the thorough study of the effects of disorder on 
strongly interacting quantum systems. Ultracold atoms 
in optical lattices m perhaps offer the most direct 
experimental system in which to realize Bose-Hubbard 
physics, however the small system sizes and destructive 
nature of many measurements limit the efficacy of exper¬ 
iments. Dimerized quantum antiferromagnets present a 
compelling alternative environment due to an exact map¬ 
ping to a lattice gas of bosons with hard-core repulsion 
EMI- These systems consist of lattices of pairs of spins 
(dimers) which, in the ground state, are all in a singlet 
configuration. This state can be viewed as an ‘empty’ 
lattice while a local triplet excitation can be thought of 
as a site occupied by a spin -1 boson (‘triplon’). 

Condensation of these bosons corresponds to exotic 
magnetically ordered states seen in materials such as 

TICUCI 3 [Him]. CS 2 CUCI 4 [I2[T3], BaCuSi 2 06 [TH1ITK] . 

SrCu 2 (B0 3 ) 2 PH 01!, and Ba 3 Mn 2 0 8 HMQj- These 
systems provide excellent experimental setups to probe 
quantum critical behavior through field- and pressure¬ 
tuning, and have motivated some notable theoretical 
works based on bond-operator techniques [ 2 T 1 BH] . 

Recent experiments on disordered quantum antifer¬ 
romagnets have seen evidence for novel glassy phases, 
particularly in bromine-doped dichloro-tetrakis-thiourea- 
nickel (DTN) J2S] where both Bose and Mott glass phases 
of bosonic quasiparticles have been observed. Such 
phases have also been seen in other materials [HHl - lHHl and 
in quantum Monte Carlo simulations [30lf34] , 

Motivated by these experimements, in this Letter we 
present an analytic treatment of dimerized quantum an¬ 
tiferromagnets with weak intra-dimer bond disorder us- 



FIG. 1: (Color online) Square lattice of dimers with intra¬ 
dimer coupling Jo. Adjacent dimers are coupled by exchanges 
Ji and J 2 . 


ing the hard-core boson formalism. We perform a strong 
coupling expansion [35) !36i combined with a replica dis¬ 
order average to derive an effective held theory. From a 
renormalization group (RG) analysis we obtain the phase 
boundaries between the gapped magnetic states - or in 
boson language, incompressible Mott insulating states - 
and the adjacent spin-glass phases. We show that away 
from the tips of the Mott lobes, the spin glass is equiv¬ 
alent to a compressible Bose glass, while at the tips we 
have strong indication for the existence of an incompress¬ 
ible Mott glass. 

The finite compressibility of the Bose glass turns out 
to be a direct consequence of replica symmetry break¬ 
ing (RSB), a mathematical property signifying the non- 
ergodic nature of the glassy states. Our work clearly 
shows that RSB in disordered Bose-Hubbard models is 
directly linked to the physical properties of the glassy 
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phases and that it finds a natural interpretation in terms 
of analogous disordered spin systems. 

We start from a Hamiltonian describing a lattice of 
coupled dimers of S = 1/2 spins, subject to single-dimer 
anisotropy D and magnetic field h, 


H 


E [ J ° 8*1 • s i2 - D{s z a + 4) 2 - m4 + 4) 


+ 2^51 Jijmn ^im ’ &jm 
izfij m,n 


(i) 


where i,j label the dimers and m,n = 1,2 the component 
spins of the dimers. This Hamiltonian is quite generic 
and describes most of the aforementioned dimer com¬ 
pounds gum HO- For simplicity, we assume that the 
dimers are located on the sites of a d-dimensional hyper- 
cubic lattice and that the couplings between dimers are 
isotropic along different bond directions. To be specific, 
we consider superexchanges J\ and J 2 between adjacent 
dimers (Fig. [I]), where Jo 3> J\ > J 2 > 0. 

The mapping of the dimerized quantum antiferromag- 
net to a model of hard-core bosons is achieved by express¬ 
ing the spin operators of each dimer in terms of singlet 
and triplet bond operators pUrfMj . s"*" | 0 ) = ^=(|tl) — 

lit)), tl |0> = |tt), 4 10 ) = ^(|tl) + lit)), and ti | 0 ) = 
|||). Using the hard-core constraint s^s + ft a t a = 1 we 
obtain 


H = - E Merrier + V E E fTl<72 

i,cr (i,j) gi ,&2 

+ 1 E [(C - - t] + ) + h.c. 

<hi) 


( 2 ) 


with h ia = i\ a iia, t = (Ji - J 2 )/2, V = (Ji + J 2 )/2, and 
Ha- = — (Jo — D) + ah. The inclusion of an anisotropy 
D > 0 is not crucial for our analysis but simplifies mat¬ 
ters by allowing us to project out the to triplet which is 
energetically unfavorable. 

The different Mott insulating states can be easily found 
in the atomic limit t —>• 0. For \h\ < Jq — D, the occu¬ 
pation numbers of both triplets are zero, corresponding 
to a gapped non-magnetic state. For sufficiently strong 
fields the magnetization is fully saturated with exactly 
one triplon on every site, m+ = 1 or to_ = 1 , depend¬ 
ing on the sign of h. We label these states as m = ±1 
Mott insulators, respectively. Between the non-magnetic 
and fully polarized states, repulsive interactions between 
triplons on neighboring sites stabilize checkerboard or¬ 
der where every second site remains empty (m = ± 1 / 2 ). 
Dimer couplings beyond nearest neighbors lead to addi¬ 
tional incommensurate states with filling fractions that 
crucially depend on the lattice geometry. In all Mott in¬ 
sulating phases, the magnetization does not change as a 
function of field, giving rise to magnetization plateaus. 

For large enough f, the system becomes a superfluid, 
corresponding to a canted XY antiferromagnet. The 



FIG. 2: (Color online) Phase diagram as a function of dimer 
coupling x and magnetic field y for 3/0 = d(Ji± J 2 )/{Jq — D) = 
1. Dashed lines show the MI/SF mean-field transitions of 
the clean system. Disorder leads to the formation of an in¬ 
compressible Bose glass (BG) between the MI and the SF, 
which turns into an incompressible Mott glass (MG) at the 
tips of the Mott lobes. Solid red lines are the MI/BG phase 
boundaries obtained from an RG analysis in d = 3 for a 
disorder strength of 5 = A/(Jo — D) = 0.3. Inset: RG 
flow of the inverse mean Io = 1/(1 + r), and relative vari¬ 
ance Ai ~ (r 2 — f 2 )/r 2 of the random-mass distribution. 
There are three fixed points: a stable MI fixed point at 
(Jo,Ai) = (0,0), a MI/SF transition at (1,0) and a MI/BG 
transition at (0, d/9). The BG/SF transition is not accessible 
in our strong-coupling approach. 


phase boundaries are obtained from a strong-coupling ex¬ 
pansion around the atomic limit. Performing a Hubbard- 
Stratonovich decoupling of the hopping term, we obtain a 
dual continuum action for the superfluid order parameter 
ip(r ,t) in space and imaginary time [55H3U] . 

S Q = f A'(k,w)|^(k,w)| 2 +u f |V>(r,r)| 4 . (3) 

J ko; J rr 

with A'(k.w) = ( k 2 — * 7 iw ± 7 2 w 2 + r) in the momen¬ 
tum and frequency domain. The mass r and the inter¬ 
action vertex u are related to the local single- and two- 
particle Green functions for the different Mott insulating 
states respectively, e.g. r = r m = 1 + 2 tdG m (uj = 0) 
m- The MI/SF mean-field phase boundaries are ob¬ 
tained by r = 0 and shown as dashed fines in Fig. [2] as a 
function of dimensionless dimer coupling x and magnetic 
field y. The phase diagram has the familiar lobe struc¬ 
ture of Bose-Hubbard models and is symmetric around 
h = 0. It has been suggested [TIHTi] that the tips of 
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the fractionally-filled lobes may exhibit first-order or su¬ 
persolid behavior, however such a question is beyond the 
reach of the present analysis. 

The coefficients of the frequency terms are given by 
derivatives of the mass coefficient with respect to the 
magnetic field, 71 = —dr/dy and 72 = — \d 2 r/dy 2 [35] . 
At the tips of the Mott lobes, the slope 71 vanishes 
and the field theory becomes relativistic, reflecting the 
particle-hole symmetry at these points. Previous works 
EM! studying other aspects of the model have largely 
neglected the frequency terms, however we retain them 
here to study how the behavior changes near the tips. 
These frequency terms will turn out to be crucial to cor¬ 
rectly describe the thermodynamics, a feature overlooked 
by the earlier studies. 

The key question we are trying to answer in the present 
work is whether the glassy phase formed in disordered 
Bose-Hubbard models is always a compressible Bose glass 
[3E] , or if the more elusive incompressible Mott glass may 
exist at the high-symmetry tips of the Mott lobes (Ml 
SB SI S3]- First seen in Id fermion systems ]37J 135] . 
the Mott glass has also been predicted to exist in the 
0(2) quantum rotor model [3HH52], which maps to the 
Bose-Hubbard model at commensurate fillings, and has 
been experimentally observed in the disordered quantum 
antiferromagnet DTN [25] . 

We focus on random mass disorder such that fi ia = 
fi a + £j. This can come from disorder in the intra-dimer 
coupling Jo, the anisotropy D or the applied field h. In 
the following, we assume that the disorder has a sym¬ 
metric box distribution of width 2A and is uncorrelated 
between different sites. For sufficiently bounded disorder 
the phase diagram retains Mott insulating regions, e.g. 
a central Mott lobe is present for 8 = A/(Jo — D) < 1. 
The Hubbard-Stratonovich transformation is performed 
in the same way, leading to disorder in all coefficients of 
the dual action. We use the replica trick [53. to obtain 
the disorder averaged free energy. The replicated action 



consists of two parts. The first contribution is simply 
n copies of the original action with disorder averaged 
coefficients 77, 77, f, and u. The second term is the 
disorder vertex, which is non-diagonal in replica space 
a, f3 and imaginary time and proportional to the variance 
of the random-mass distribution g = (r 2 — f 2 ). 

To determine the phase diagram in the presence of 
weak disorder, we use a momentum-shell RG approach. 
As in previous work | |38H40] , we make the change of vari¬ 
ables Iq = 1/(1 +r), and introduce the relative disor¬ 
der variance A = I^g to distinguish between the Mott 
insulating and glassy phases. In all but the superfluid 
phase, Iq flows to zero, reflecting the short-ranged super¬ 
fluid correlations. The relative variance A compares the 


shift of the random mass distribution with its spread. In 
the Mott insulator, the distribution shifts faster than it 
spreads and A renormalizes to zero. If the spread is faster 
than the shift, the tail of the distribution pushes through 
zero, indicating a glassy phase where the physics is dom¬ 
inated by rare superfluid regions. Taking into account 
the 1-step RSB in this model m , the RG equations are: 


/'(f) = (3/2A 1 ^2)/ 0 + 2/ 2 , (5a) 

Ai(f) = (4/ 0 - d)Ar + 9A?, (5b) 

= (2 — z + Ai)77, (5c) 

72 CO = (2 — 2z + Ai)y 2 + Ai/q 7 i, (5d) 


where z is the dynamical critical exponent and Ao = 0 
and Ai = I 2 gi denote the step heights of the Parisi dis¬ 
order function muss]- We neglect the u vertex since it 
is irrelevant away from the tips (77 7^ 0). The RG flow in 
the /0-A1 plane is shown in the inset in Fig. [2] The MI/SF 
fixed point is unstable against disorder, confirming that 
even for infinitesimal disorder, the transition from the 
Mott insulator is into a disordered insulating state and 
not into a superfluid mi Eg. As we show later, the dis¬ 
ordered state is a compressible Bose glass, except for the 
tips where we see strong indications for an incompressible 
Mott glass state. 

The relative variance Ai diverges in the Bose glass. We 
stop integration at a scale f* where Ai(f) = 1 and our 
RG becomes invalid. This scale can be identified with 
a correlation length £ ~ which corresponds to the 
typical separation of superfluid regions. Note that £ is 
not the superfluid correlation length, which remains finite 
at the MI/BG transition. Linearizing near Pmi/bGj we 
find the correlation length diverges as £ ~ (a; — a; c ) _1 / d . 

From the dependence of Iq and Ai on the microscopic 
parameters, we can determine the MI/BG phase bound¬ 
ary as a function of the dimensionless dimer coupling 
x and magnetic field y for a given disorder strength <5 
(Fig. §- As in the conventional Bose-Hubbard model, 
the Mott insulating regions shrink when disorder is 
added, and a glassy phase intervenes between the Mott 
insulating and the superfluid regions. In the original 
magnetic language, this corresponds to a spin glass phase. 

To determine the nature of the glassy phase, we cal¬ 
culate the compressibility. For a single component bo¬ 
son system, the local compressibility is defined as the 
derivative of the local density with respect to the chem¬ 
ical potential and related to the particle-number fluctu¬ 
ations, Ki = d(n.i)/dfi = PJ2j{{nihj) — (fij)(hj )). This 
is easily generalized to the two-triplon case, where the 
compressibility is defined as the derivative of the local 
magnetization with respect to magnetic field. Perform¬ 
ing a disorder average we can calculate k = —d 2 jF/dy 2 
directly from the replica field theory, yielding ]3T l 


2 S d I 2 j 2 X 1 

{2ir) d (J 2 7 2 + 4/o7 2 ) 3 / 2 ’ 


( 6 ) 
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FIG. 3: (Color online) (a) Power-law behavior of the com¬ 
pressibility k(1*) ~ C(x—x c ) 1 2 (d = 3) of the Bose glass close 
to the transition to the Mott insulator. The blue dashed lines 
are guides to the eye with gradient 1/2. Near y = 0 there is 
an anomalous suppression of the range of universal behavior, 
(b) C versus y, showing the vanishing compressibility close to 
the tip of the central lobe at y = 0. 


where Sd is the surface of a d dimensional unit sphere. 
It is crucial here to include the effects of 1-step RSB, 
otherwise n vanishes identically as a consequence of the 
frequency structure that is inherent to all Bose-Hubbard 
models mi mg. This remarkable result reveals the 
surprising physical importance of RSB in this system. 
Linearising around Fmi/bg we find that at the transition 
to the Mott insulator the compressibility vanishes as 

k = k(P) = C(x — x c ) a_s , (7) 

where C ~ |'y 1 1 - Away from the tip, the disordered state 
adjacent to the Mott insulator is a compressible Bose 
glass. Approaching a particle-hole symmetric Mott lobe 
tip, the coefficient C vanishes, suggesting a change of 
universality. The analytical result ([7| is confirmed by a 
numerical calculation of k(£*) from the full RG equations 
© for the central Mott lobe in d = 3 with S = 0.3 (see 
Fig. |3}i). The coefficient C has a strong dependence on 
the field and vanishes linearly as y —>■ 0 (Fig. ©>)• Very 


close to the tip at y = 0, there is a strong suppression 
of the range of universal behavior, but no indication of a 
crossover to a different universality class. This highlights 
the singular nature of this point and strongly indicates 
the existence of an incompressible Mott glass state. 

Although we have focused on dimerized quantum an- 
tiferromagnets using the hard-core boson formalism, our 
prediction of a Mott glass is valid across a wide range 
of systems, including conventional Bose-Hubbard and 
Jaynes-Cummings Hubbard models. Any such Bose- 
Hubbard-like model may be treated using the methods 
outlined here, with the microscopic differences appearing 
only in the UV-scale starting values of the flow parame¬ 
ters. This prediction lends weight to previous numerical 
quantum Monte-Carlo work mi ezh. It may also ex¬ 
plain the controversy over the existence of a direct MI/SF 
transition at the tips of the Mott lobes in the disordered 
Bose-Hubbard model: previous works which used com¬ 
pressibility as the criterion for the onset of a glassy phase 
will necessarily have missed the transition between the 
Mott insulator and the Mott glass. 

The breakdown of self-averaging [351 EH E2] , the im¬ 
portance of replica symmetry breaking [40] , and the con¬ 
nection with spin-glass phenomena strongly suggest that 
the formation of Bose and Mott glass phases is not sim¬ 
ply a weak localization phenomenon. It would be inter¬ 
esting to review the nature of these phases in the context 
of many-body localization ,59] and the related entangle¬ 
ment entropy scaling [60 . 

The equivalence between Bose-Hubbard models and 
dimerized quantum antiferromagnets allows for multiple 
complementary experiments to verify our theoretical pre¬ 
dictions. Previous measurements on disordered ultracold 
atomic gases have inferred the presence of the Bose glass 
from macroscopic measurements |61K63| , however the bo¬ 
son number fluctuations associated with the local com¬ 
pressibility are also now within reach of quantum gas 
microscope systems SUE potentially allowing for direct 
imaging of the glassy phases. Thermodynamic measure¬ 
ments of the magnetic susceptibility and the specific heat 
of dimerized quantum antiferromagnets can provide clear 
signatures of Bose and Mott glass phases £25]. We also ex¬ 
pect to see characteristic differences in the glassy dynam¬ 
ics, which could be studied with /rSR, as well as in the 
magnetic excitation spectra. Dimer systems exhibiting 
geometric frustration are a particularly intriguing theo¬ 
retical problem for further study, as are additional types 
of disorder such as non-magnetic impurities 133] , 
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Supplementary Material 


auxiliary Hubbard-Stratonovich fields ipi(r), 4 > i ( T ), 


In this Supplementary Material to our Letter we go 
through the technical details of our calculation. We show 
in detail the construction of the strong coupling field the¬ 
ory and the calculation of the compressibility. In particu¬ 
lar, we wish to stress the importance of including replica 
symmetry breaking effects to correctly describe the dis¬ 
ordered phases. 


Mapping to Hard-Core Bosons 


The spin operators of each dimer may be written in 
terms of the singlet and triplet bond operators as: 


qo. _ 

*i,2- 2 


•f'.sVrt ± /?,.s it a (3 'ytpt'-) 


( 8 ) 


where |o> = ^(|tl) - lit)), 4 |0) = -^(|tt) - III)), 

*110) = + III)), and P z |0) = ^(|tl) + lit))- 

These operators are bosonic and subject to the hard-core 
constraint s' s + i' a i a = 1- The local dimer Hamiltonian 
is diagonal in the standard triplet basis, t\_ |0) = |tt) = 
|1,1), H |0) = III) = |1,-1), and 4 |0) = t\ |0) - |1,0). 
The Hamiltonian (Eq. 2) in our Letter is obtained by 
a change of variables to this basis. Note that we have 
also used the hard-core constraint to eliminate the sin¬ 
glet operators and to rewrite on-site products of number 
operators, = h and fii_hi + = 0. 


Strong Coupling Field Theory 

We begin from the bosonic Hamiltonian shown in our 
Letter (Eq. 3) and use it to write the partition function 
at T = 0 as an imaginary time path integral over bosonic 
coherent states, Z = f T>[(p,(p]e~^ So+St \ where the local 
and kinetic contributions to the action are given by 


e -S t [4>M 
S t [ip, ip] 

c' 


J IMl+ s ;»>, (10a) 

/£ (10b) 



[0»_(t) - ^+(r)]V>i(r) + h.c. (10c) 


Because of the structure of the initial action, we only 
require a single Hubbard-Stratonovich field to perform 
the decoupling, despite the presence of two bosonic fields 
in the initial action. We can now trace out the original 
fields to obtain the partition function in terms of the new 
fields, Z = J T>[ip,ip]e~^’^h The resulting dual action 
is given by S[ip,ip\ = St\ip , ip] — SS[ip, ip] with 


6S = In /Texp j J + 6|(r)^(r)]| \ 11) 

where T denotes the time-ordering operator and the av¬ 
erage is taken with respect to the local Hamiltonian, 

W-0 ^ ) f-lall'ia T H y ' y [ CJi<T2nj ( 7 1 ?tj ( 7 2 , (12) 

i,cr=± (i,j) a h a 2 


flier = For brevity, we have defined 4(r) = 

ii- (t) - 4 + (r) and f icr (r) = e HoT t i(T e~' HoT . In principle, 
we can expand 5S to any desired order in the fields ip. 
However, since Wo contains interactions, we cannot use 
Wick’s theorem to calculate expectation values of prod¬ 
ucts of the original bosonic operators. The expansion 
only contains powers of \ip\ 2 since Wo preserves the local 
particle number. Expanding SS up to quartic order, we 
obtain 


^ = E / GMI^MI 2 + E / l^( T )l 4 ’ ( 13 ) 

J UJ j T 


where G(u>) is the Fourier transform of the single particle 
Green function 


_ r°° c 

So [<M] = / dr | E (^i<r( T )(^ - Vo)4>i<j{T) 

J o 1 i,*=± 

+v E E CTifT2 l^lCTl i T ) | I0JIT2 (r)| 2 },(9a) 

(i,j) <7 ii G 2 

nOO 

St [</>, 4>]= dr E T ij \4>i -( r ) - <fo+ ( T )] 

x [<t>j-( T )~4>j+( T )\- ( 9b ) 

The hopping matrix elements Xj, are equal to t if i,j 
are nearest-neighbor sites and zero otherwise. To facili¬ 
tate the strong-coupling expansion around the local limit 
(t —> 0), we decouple the kinetic term by introducing 


G(u) = - dre < " T (TSi(0)SJ(r))o. 


(14) 


It is straightforward to calculate G m (uj) for the different 
MI states, 


G±i(lj) = - 
^±1/2 (^) = — 


1 


li± — zV + ico ’ 

1 


1 


+ iu> — n± + zV — iio 

1 


GoM = 


Hzf + zV — iuj' 

1 1 


fi+ + iui i-i- — iu 


(15a) 

(15b) 

(15c) 



The coefficient U of the quartic vertex is given by the 
connected parts of the two-particle Green functions. Un¬ 
like in the conventional Bose-Hubbard model, the calcu¬ 
lation here is relatively simple because of the hard-core 
constraint. In each of the regions, U m is given by 


U ± i 


u± 1/2 


U 0 


1 


(/G - 2 dVY 


2 {fx± — 2 dV) 2 (g± — /Zip — 4 dV) ’ ^ ^ 

11 1 
/4 + {2dV-^ ± f ~ (, n T + 2d.Vf 
AG + AG 


(fj,± — 2dV) 2 (n^ + 2dV) 2 
1 1 
2/4(AG “AG)’ 

1 1 -I- 

a4 


(16b) 

(16c) 


Analyzing the above expression we find that U is al¬ 
ways positive. Hence there is no indication of a first-order 
instability of the m = ± 1/2 checkerboard lobes at this 
level. The hopping term St can be treated in the long- 
wavelength limit. In this case it is possible to invert the 
original hopping matrix and we obtain 


S t 




IV’(k, W )| 2 . 


(17) 


Combining St and the quadratic contribution to 5S 
and expanding G(w) for small frequencies we obtain 


S 2 = J + R ~ iK i u + K 2 u 2 ^J |^(k, w )| 2 

(18) 

with R = 1/2 dt + G(u> = 0). Inspecting G(u>) and re¬ 
membering that = —(Jo — D) + ah, we see that G is 
a function of h + iui. From this, it immediately follows 
that Ki = —dR/dh and K 2 = — \d 2 R/dh 2 . 

As a last step, we rescale the action to dimension¬ 
less units, using w —A Vw, k —► V2dk, and \ip\ 2 —A 
y/2d t/V-\ij)\ 2 . This leads to the dimensionless strong¬ 
coupling action given in the Letter with coefficients 
r = 2 dtR, = 2 dtVK 1 , q 2 = 2dtV 2 K 2 , and u = 

V2t~ 3d t 2 U/V 3 . 

In terms of the dimensionless dimer coupling (hopping) 
x = d(J\ — J 2 )/(J 0 — D) and magnetic held (chemical 
potential) y = h/(Jo — D) the mass coefficients on the 
different regions are 


ro = 1 — a:[(l + y) 1 + (1 — 2/) _1 ], (19a) 

r± 1/2 = 1 — ^[(1 + UoT y) 1 — (1T 2/) _1 

+{l~y 0 ±y)~ 1 }, (19b) 

r ±i = 1 + x(l + y 0 ± 2/) _1 , (19c) 


Replica disorder average 


In the presence of on-site disorder e*, which randomly 
shifts the local chemical potentials, + et, the 

Hubbard-Stratonovich transformation can be performed 
in exactly the same way. This results in a dual strong¬ 
coupling action on a lattice with disorder in all the coef¬ 
ficients, e.g. the mass on site i is given by 

n = r{n a + €i,V,t), (20) 

where r(p a , V, t) is the the mass expression we have de¬ 
rived for the clean, homogeneous system. As discussed 
in the Letter, we use the replica trick J- = —Tin Z = 
—Timing \ — 1) to average the free energy over the 

disorder. Here we provide some additional details. The 
effective replica action <S e ff is obtained by taking n copies 
of the system (labelled by a) and performing the disorder 
average 

Z™ = J V\$ a ,'4) a \e-'Z<* s ®on'l’o.} 

=: j V^ a ,i, a }e- S ^U’M\ (21) 

Expanding the exponential, performing the disorder 
average, and re-exponentiating, we obtain the general ex¬ 
pression 

5 c ff = ^2 S a — - ^2 (SaS/3 - S a S/3 ) , (22) 

at otfi 


where S a = S[ip a ,ip a ], for brevity. This expression is 
exact up to quartic order in the fields. The first term 
simply corresponds to n copies with disorder-averaged 
coffiecients r, 7 1; y 2 , u. For a box distribution p(e) all 
the averages can be performed analytically. E.g, for the 
central Mott lobe around zero field we obtain the average 
mass 


x m +s) 2 -y 2 -\ 
25 [(1 — <5) 2 - y 2 _ 


(23) 


with x the dimensionless hopping, y the dimensionless 
field, and 5 the dimensionless width of the box distri¬ 
bution, as defined in the Letter. Since the field deriva¬ 
tives commute with the disorder average, the relations 
7 y 1 = — df/dy and 7 2 = —\d 2 r/dy 2 still hold. It is im¬ 
portant to include the moments of p(e) to all orders in 
order to account for the boundedness of the distribution. 
This is crucial for the presence of Mott insulating states 
in the disordered system. 

The second term in Eq. (22) contains various quartic 
disorder vertices. We neglect all the vertices with tem¬ 
poral gradients as they are irrelevant under the RG, and 
keep only the disorder vertex that arrises from the repli¬ 
cation of the mass term. It is obvious from Eq. (221 that 
the coefficient is given by the variance g = r 2 — f 2 of the 
random mass distribution. This can also be calculated 
analytically for the box distribution p(e). 


where y 0 = d(Ji + J 2 )/(Jo - D). 
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Renormalization and Replica Symmetry Breaking 

We integrate out the highest energy modes, corre¬ 
sponding to momenta in the infinitesimal shell e~ dl < 
|k| < 1. To regain resolution, we rescale momenta 
k —>■ ke _d/ and frequencies u —> uje zdl . The fields are 
rescaled such that the k 2 term remains constant. This 
leads to coupled RG equations for r, qq, 7 2 , u , and g. 
For simplicity, we neglect the u vertex since it is irrel¬ 
evant under the RG and only leads to a small shift of 
the phase boundaries. We allow for replica symmetry 
breaking (RSB) and replace g by a matrix g a p of general 
Parisi hierarchical form |53j . In previous work uni we 
examined this model and showed the presence of 1 -step 
RSB. Here we will show that RSB is responsible for the 
finite compressibility of the Bose glass. 

Taking the replica limit, the matrix g a p becomes a 
monotonically increasing function g{u) (u £ [ 0 , 1 ]), which 
describes how the off-diagonal matrix elements of g a p de¬ 
crease as we move away from the diagonal with elements 
9aa = 9- As we have discussed in the Letter, the relative 
variance is the relevant variable to distinguish the Mott 
insulator and the Bose glass. We therefore rescale the 
Parisi matrix by Iq, [A,A(it)] = Iq [g, g(u)]. Note that 
I 0 = 1/(1 + r) is the on-shell propagator in the zero- 
frequency limit and for large r (deep in the insulating 
phases) becomes the inverse mean of the random mass 
distribution. At one-loop order the RG equations are 
given by m 


solutions have 1-step RSB. To be more precise, A (u) 
is a step function with A (it < u c ) = Ao = 0 and 
A (it > u c ) = Ai = A where the step position is u c = l/2po 
HO] . Specifying to this step function, we obtain the set 
of RG equations given in the Letter. 


Compressibility 

In this section we will provide details of the compress¬ 
ibility calculation and show that the finite compressibility 
of the Bose glass is a consequence of RSB. The compress¬ 
ibility is defined as 


d 2 T 1 d 2 Z n 

~W ~ n~d^~' 


(26) 


where y denotes the dimensionless magnetic field. Tak¬ 
ing into account that all the coefficients in the quadratic 
replica action are functions of y we obtain 


= lim n -( 2 7 f (IVg^aw)! 2 ) 
n ^ 0n y 

+ ^2 [ (7i + 2 *^7 2 ) (7i + 2 iw , 7 2 ) 

a p J kk'ww' 

x (|V. Q (fc, W )| 2 |^(fc', W ')| 2 ) 


(27) 


d/p 

d£ 

d7i 

d£ 

dj 2 

d£ 

dA 

d£ 


dA(it) 

d£ 


—2 + A + po ^A — (A)^ 

Io + 2I 2 0 , 

(24a) 

(2 — z + A) 7 i, 


(24b) 

2 — 2 z + A)y 2 + A/q 7 i, 


(24c) 


(4/p -d) A + 6 A 2 


+2po ^A — (X)J A + 2 ^A 2 — (A 2 )^ 


(4 1 0 — d)X(u) + 2A 2 (u) + 4 AA(m) + 4p 0 


Na- 


(24d) 


A(it) — dv [A(m) — A(i>)]^ 

Jo 

(24e) 


where the averages are taken with respect to the full ac¬ 
tion with interaction and disorder vertices u and g a p. 
Note that we have kept terms up to quadratic order 
in frequency as we did in the derivation of the strong 
coupling action. We compute the expectation values 
by expanding to linear order in the vertices, (...) ~ 
(.. .)o — (... S u )o — (... S g ) o, which results in three addi¬ 
tive contributions to the compressibility, k = kq+K u +K g . 
Note that the replica limit n —► 0 ensures that only con¬ 
nected diagrams contribute. 

Let us first calculate the compressibility kq for u = 
gap = 0. Taking the expectation values with respect 
to the quadratic action, (...) = (...) 0 , we immediately 
obtain 


K 0 = 


_ _ 2 d _ d o ® 

2 % - 7 ,^+ 4717 ,^ + 47 , ^ 


h, (28) 


where (A 71 ) = X n (u)dn and 

Po = 4 _1 ( l + r-*7iw + 7 2 w 2 ) -1 

J UJ 

1 

\/ Jill + 4J o7 2 


where I\ denotes the integral of the correlation 
function over momentum and frequency, I\ = 
fkcj (^ 2 — irf lU + 7 2 ^ 2 ) • Since in the insulating 

phases f —> oo under the RG, we can approximate the 
momentum integral for |k| < 1 <C r, and carry out the 
frequency integration to obtain 


Analysis of the RG equations (24) shows that the sys¬ 
tem is unstable towards RSB and that the only stable 


S d 1 

(2n) d -|- 4(1 -|- f)^ 2 


( 29 ) 
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with Sd the surface of the d-dimensional unit sphere. 
Taking the derivatives in Eq. (28), we find that kq van¬ 
ishes identically. The same is true for k u since it is pro¬ 
portional to further derivatives of k 0 . This simply reflects 
that the gapped MI is incompressible. 

A finite compressibility has to be a consequence of the 
diverging disorder vertex, n = n g . Starting from a gen¬ 
eral Parisi matrix g a p, we obtain 


\ VIlil + hi 2 the step position. Substituting into 
Eq. ( |30| we arrive at the final expression provided in 
our Letter. 


Critical Scaling 


s d 

(2ny 


A- 


4/ 0 2 7l 


( 7 o7? + 4o7 2 ) 2 ’ 


(30) 


where [A, A(u)] represents A a p = I^ga /3 after taking the 
replica limit, as we have explained above. This result is 
remarkable. It demonstrates that without RSB, g a p = g 
or A (u) = A, the compressibility vanishes. This demon¬ 
strates that RSB is essential to correctly describe the 
physical properties of the Bose glass phase. 

For the stable 1-step RSB solution, we obtain A — 
(A) = u c Ai with Ai the step hight and u c = l/2po = 


Using a scale-dependent dynamical exponent z(£) = 
2 + Ai(£), we enforce that ^ does not scale under 
the RG. In the vicinity of Pmi/bg we obtain 7 2 (£) = 
|c 0 7 i exp[— (2 — d/6)£]+Ci exp[— (2+d/9)£]. If we are not 
exactly at the tip (qq 0), the first term will dominate 
for sufficiently large i and q 2 (0 ~ 7 i/o(Z). Substituting 
into Eq. (6) in the main text and evaluating at the cor¬ 
relation length £ ~ e l ~ (x — x c )~ 1 ^ d , we find that at 
the transition to the Mott insulator the compressibility 
vanishes as k = k(£*) = C(x — x c )a _ B as stated in our 
Letter. 






